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SUMMARY 


Stability in energy for the Newmark 3-family of time integra- 
tion operators for nonlinear material problems is examined. It 
is shown that the necessary and sufficient conditions for uncondi- 
tional stability are equivalent to those predicted by Fourier 
methods for linear problems. 


INTRODUCTION 

In this paper, stability in energy for the Newmark family 
(ref. 1) of time integration operators is examined. Stability for 
these operators was considered in the original paper of Newmark, 
who used essentially Fourier techniques which are strictly appli- 
cable only to linear problems. Belytschko and Schoeberle (ref. 2) 
have shown the unconditional stability of the particular form of 
the Newmark 3 - operator that corresponds to the trapezoidal rule 
(y=%, 3=%) for nonlinear material problems by energy methods. 

Hughes (ref. 3) extended this proof to the range of parameters 
( y—h , 3-%) • In this paper, it is shown by generalizing the defini- 
tion of discrete energy, sufficient conditions for unconditional 
stability in energy on both y and 6 can be obtained. These condi- 
tions are equivalent to the necessary conditions for the uncondi- 
tional stability of the Newmark operators in linear problems, so 
the conditions obtained herein are necessary and sufficient for the 
unconditional stability for nonlinear material problems. 


PRELIMINARY EQUATIONS 

The equations will here be presented in the formalism of the 
finite element method. As indicated in Belytschko, et al (ref. 4), 
the spatial discretization in finite difference methods is ba- 
sically identical, so the choice of finite element notation .is 
only a matter of convenience, not a restriction on the proof. The 
equations will only be outlined; details may be found in Zienkiew- 
icz (ref. 5) . 

The fundamental step in any spatial discretization, which is 
often called the semidiscretization/ is a separation of variables 
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in the form 


u(x,t) = £(x) (t) (1) 

where x is the Cartesian coordinate, t the time, u the displace- 
ment fi&ld, £ the shape functions, and d the nodal displacement 
of element e. The strains can then be related to the nodal dis- 
placements by 

e * Bd (e) = B L (e) d (2) 

( 0 ) , 

where B consists of derivatives of the shape functions and L ' is 
the connectivity matrix. The discrete equations of motion are 
then 


M a + f = £ 


(3) 


where M is the mass matrix, a the nodal accelerations (second de- 
rivatives of d with respect to time) , £ the external nodal forces 
and f the internal nodal forces, which are given by 

1=2 L (e)T f (e) = 2 i (e)T / § T o dv (4) 
e e v( e ^ 


Equations (3) and (4) can be derived from the principle of virtual 
work with the inertial forces included in a d'Alembert sense; see 
for example Belytschko, et al (ref. 6) . 

We define a discrete internal energy by 



where upper case subscripts denote the time step and A denotes a 
forward difference 


ASj = £ I+1 - (6a) 

and 

0 < ij < 1 (6b) 
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When y=%, eq. (5) represents a trapezoidal integration of the non- 
linear stress strain curve, while y=0 corresponds to Euler integra- 
tion. 

By means of eqs. (2) and (4) , eq. (5) can also be written as 

u i+i = u i + 35 Adj [d-yjfj + yf 1+1 ] (7) 

We require that the discrete internal energy be positive definite, 
so that 


I t 

£ AejCd-yJdj + y£j+i 1 - 0 

J=1 


The kinetic energy T is given by 

T 

T I = 35 -I - -I 


(9) 


where v are the nodal velocities (first derivative of d with res- 
pect to time) . 

The Newmark difference formulas are 

v I+i = y_x + At [(l-yJaj + ya I+1 ] (10) 

^X+l = dj + At Vj + At 2 [(%-3)a I + 3a I+1 ] (11) 

When 3 > 0, these equations are implicit, and hence for nonlinear 
materials, the solution of a nonlinear system of equation is 
necessary. The exact solution of the nonlinear system of equations 
at each time step is not possible; at each time step there will be 
an error f| rrl given by 

f j rr = p x - f x - Ma r (12) 

We define an energy error criterion 

|Ad£ [(l-y) f| rr + < e (u i + T i ) (13) 

where e is a small constant and require that the solution of the 
nonlinear equations at each time step satisfy this criterion. 
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PROOF OF UNCONDITIONAL STABILITY 

We will now show that the error criterion, eq. (13) , is a 
sufficient condition for the unconditional stability in energy of 
the Newmark integration formulas for y^h, 3 - Y/2. Stability in 
energy is described in Rich tiny er and Morton (ref. 7) and has pre- 
viously been used for the derivation of stability conditions for 
the solution of linear problems by the Newmark 0-method by Fujii 
(ref. 8) . 

The demonstration of unconditional stability in energy requires 
that it be shown that a positive definite norm of the solution is 
bounded regardless of the size of the time step. As pointed out 
in reference 7, the norm need not be the physical energy, though 
in many cases it is. For the purposes of this proof, we define 
the norm by 

S z = T- x + Uj + (20-Y)a I T Ma I (14) 

Because of the requirement of eq. (8) , Uj is positive definite, 
whereas the positive definiteness of the mass matrix M assures the 
positiveness definiteness of the remaining two quantities if 


20 - y (15) 

Stability in energy is then assured if we can show that Sj is 
always bounded, i.e. that 

S I+1 < (l+e*)S I (16) 


where e* is an arbitrarily small quantity. The interpretation of 
the condition of eqs. (14) and (16) is as follows. Provided that 
the discrete internal energy is a monotonically increasing func- 
tion of the displacements, the boundedness of Si implies that the 
velocities and displacements are bounded, which corresponds to the 
notion of stability. 

The proof of stability then consists of deducing eq. (16) from eqs. 
(7) to (11) and the homogeneous form of eq. (3) . From eqs. (9) and (10), 
it follows that 

T I+1 = T I + At If 1 -*)®! + Ya I+1 ] 


At' 


[d-YJa/ + ya^] M [(1-Y)a.j + ya I+1 ] 


T 

Y -I+l 


(17) 
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so 


T 


1+1 


Tj. + [Adj + At 2 ($-%Y)a I T + At 2 (%y-6) a^-J 

«[a-Y)» x + rs I+1 ] 


Thus if we let 


(18) 


Y = U 

then eqs . (7) and (18) yield 

T I+1 + U I + 1 = T I + u i + Ad£[(l-Y) (Ma I+ f T ) 

+ y(Ma I+ ^ + fj + ^) J + At “ 2 )— i + '('2 I+l] — 

[(l-y)^ + ya I+1 ] 


(19) 


( 20 ) 


The second term on the right hand side of eq. (20) corresponds to 
the error in energy as defined by eq. (13) , and the last term can 
be rearranged so that we obtain 

T I+ i + U I+1 1 (l+eJlTj+Uj) +4^7-23) (a^Ma^ - ^ M a J .) 

+ At 2 (y-3s) (y- 28) (a^ x - a*) M (a I+1 - (21) 

The last term on the right hand side of eq. (21) is negative semi- 
definite if 


Y > h and 28 £ y (22) 

or if both of the above inequalities are reversed. However, if the 
inequalities are reversed, as can be seen from eqs. (14) and (15) , 
the norm Sj is not positive definite. Hence, only the conditions 
given by eq. (22) are pertinent. Under these conditions, the in- 
equality of eq. (21) applies even if the last term is dropped. The 
remaining terms then yield eq. (16). 
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DISCUSSION AND CONCLUSIONS 


Several remarks should be noted in applying these results to 
computations. First, the stability hinges on the achievement of a 
solution in each time step that satisfies eq. (13) . The conver- 
gence of solution schemes, such as the modified Newton Raphson 
method, cannot be assured, and is therefore the primary obstacle 
in obtaining stable solutions. . The difficulties are particularly 
severe in elastic-plastic problems if the tangential stiffness 
method is used whenever unloading takes place ovet a large part of 
the mesh. 

It is also not clear whether the form of the error criterion, 
eq. (13) , is suitable for very fine meshes. Numerical experiments 
indicate that it becomes increasingly difficult to satisy eq. (13) 
for finer meshes, for although the criterion appears to be mesh- 
independent in that the right hand side increases with the size of 
the physical problem, the right hand side does not vary as a mesh 
is refined. Furthermore, in very large meshes there is a possi- 
bility of cancellation of errors, i.e. positive error energy 
transfer in one portion, with negative error energy transfer in 
another portion. This can be avoided by placing the absolute 
value within the summation. 

Results have been reported for a special case of this operator 
(y=h , $=%) in reference 2. Both material and geometric nonlinear- 
ities were included in those problems. However, the proofs given 
here and in reference 2 require the absence of geometric nonlin- 
earities; if geometric nonlinearities are included, eq. (5) does 
not imply eq. (7) , for in geometrically nonlinear problems AB does 
not vanish. Hence, as shown in reference 9, in geometrically non- 
linear problems, energy transfer is associated with the rotation 
of a stressed member: this effect results in the generation of 
energy if the stress is tensile and is hence destabilizing under 
those conditions. In many structural dynamics problems, the total 
rigid body rotation that takes place is insufficient for this 
energy generation to be significant. However, test problems have 
been devised where the energy error is so large that for practical 
purposes the computation can be considered unstable. 

Finally, we comment on some experience with the requirement 
of eq. (8) . This condition requires that the numerical integra- 
tion of the internal work always yield a positive quantity. In 
elastic-plastic materials and other strongly dissipative materials, 
this condition poses no problems. However, when the stress is a 
single-valued function of the strain, eq. (8) can easily be 
violated in cyclic load paths. However, numerical experiments do 
not indicate that violation of eq. (8) results in any catastrophic 
failure of the computation. 
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